#include "index_gpu.cuh"
#include "boundary_gpu.cuh"
#include "Global_Variables_extern_gpu.cuh"
// #include "Global_Variables_gpu.cuh"

//================================================================================================================================================================= =
//----------------------off-wall bounce back for wall boundary----------------------
// currently, only the domain boundary is no-slip wall, there is no particle in the fluid domain
// for one-step index (OSI) algorithm
//================================================================================================================================================================= =
__global__ void wall_BB_BC_GPU_face_xy(I_INT imax, T_P *pdf_old_gpu, I_INT pitch, int ntime)
{
    int i_f; // column //index of discrete velocities
    // int i_grid ;    //row  //index of grid position coordinates
    T_P *p_f; // the point to rowHead(fi in first grid)

    // Indexing (Thread)
    I_INT i = threadIdx.x + blockIdx.x * blockDim.x;
    I_INT ix, iy, iz;
    I_INT i_Global; // the index of point in the global distribution function array
    // Represents the position in the grid that needs to be solved
    // int position;
    T_P f[NDIR];
    if (i < imax)
    {
        // z direction (top wall)
        ix = i % (NXG0_D) + 1;
        iy = i / (NXG0_D) + 1;
        iz = NZG0_D;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = BB_xy_top_gpu[i3];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            f[i_f] = p_f[p_step_index_D((i_Global + Length_gpu[i_f]))];
            p_f = (T_P *)((char *)pdf_old_gpu + Reverse_gpu[i_f] * pitch);
            p_f[p_step_index_D(i_Global)] = f[i_f] - prc(6.) * w_equ_gpu[i_f] * Density_gpu * Lid_velocity_gpu * (T_P)ey_gpu[i_f]; // - prc(6.) * Density_gpu * w_equ_gpu[i_f] * prc(0.5)*dt_LBM_gpu*body_accelerate_gpu[2] * (T_P)ez_gpu[i_f];
        }
        //__syncthreads();
        // z direction (bottom wall)
        iz = 1;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = Reverse_gpu[BB_xy_top_gpu[i3]];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            f[i_f] = p_f[p_step_index_D((i_Global + Length_gpu[i_f]))];
            p_f = (T_P *)((char *)pdf_old_gpu + Reverse_gpu[i_f] * pitch);
            p_f[p_step_index_D(i_Global)] = f[i_f]; //- prc(6.) * Density_gpu * w_equ_gpu[i_f] * prc(0.5)*dt_LBM_gpu*body_accelerate_gpu[2] * (T_P)ez_gpu[i_f];
        }
    }
}
__global__ void wall_BB_BC_GPU_face_yz(I_INT imax, T_P *pdf_old_gpu, I_INT pitch, int ntime)
{

    int i_f; // column //index of discrete velocities
    // int i_grid ;    //row  //index of grid position coordinates
    T_P *p_f; // the point to rowHead(fi in first grid)

    // Indexing (Thread)
    I_INT i = threadIdx.x + blockIdx.x * blockDim.x;
    I_INT ix, iy, iz;
    I_INT i_Global; // the index of point in the global distribution function array
    // Represents the position in the grid that needs to be solved
    // int position;
    T_P f[NDIR];
    if (i < imax)
    {
        // x direction (front wall)
        iy = i % (NYG0_D) + 1;
        iz = i / (NYG0_D) + 1;
        ix = NXG0_D;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = BB_yz_front_gpu[i3];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            f[i_f] = p_f[p_step_index_D((i_Global + Length_gpu[i_f]))];
            p_f = (T_P *)((char *)pdf_old_gpu + Reverse_gpu[i_f] * pitch);
            p_f[p_step_index_D(i_Global)] = f[i_f]; //- prc(6.) * Density_gpu * w_equ_gpu[i_f] * prc(0.5)*dt_LBM_gpu*body_accelerate_gpu[2] * (T_P)ez_gpu[i_f];
        }
        // x direction (back wall)
        ix = 1;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = Reverse_gpu[BB_yz_front_gpu[i3]];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            f[i_f] = p_f[p_step_index_D((i_Global + Length_gpu[i_f]))];
            p_f = (T_P *)((char *)pdf_old_gpu + Reverse_gpu[i_f] * pitch);
            p_f[p_step_index_D(i_Global)] = f[i_f]; //- prc(6.) * Density_gpu * w_equ_gpu[i_f] * prc(0.5)*dt_LBM_gpu*body_accelerate_gpu[2] * (T_P)ez_gpu[i_f];
        }
    }
}
__global__ void wall_BB_BC_GPU_face_xz(I_INT imax, T_P *pdf_old_gpu, I_INT pitch, int ntime)
{

    int i_f; // column //index of discrete velocities
    // int i_grid ;    //row  //index of grid position coordinates
    T_P *p_f; // the point to rowHead(fi in first grid)

    // Indexing (Thread)
    I_INT i = threadIdx.x + blockIdx.x * blockDim.x;
    I_INT ix, iy, iz;
    I_INT i_Global; // the index of point in the global distribution function array
    // Represents the position in the grid that needs to be solved
    // int position;
    T_P f[NDIR];
    if (i < imax)
    {
        // y direction (right wall)
        ix = i % (NXG0_D) + 1;
        iz = i / (NXG0_D) + 1;
        iy = NYG0_D;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = BB_xz_Right_gpu[i3];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            f[i_f] = p_f[p_step_index_D((i_Global + Length_gpu[i_f]))];
            p_f = (T_P *)((char *)pdf_old_gpu + Reverse_gpu[i_f] * pitch);
            p_f[p_step_index_D(i_Global)] = f[i_f]; //- prc(6.) * Density_gpu * w_equ_gpu[i_f] * prc(0.5)*dt_LBM_gpu*body_accelerate_gpu[2] * (T_P)ez_gpu[i_f];
        }
        // y direction (left wall)
        iy = 1;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = Reverse_gpu[BB_xz_Right_gpu[i3]];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            f[i_f] = p_f[p_step_index_D((i_Global + Length_gpu[i_f]))];
            p_f = (T_P *)((char *)pdf_old_gpu + Reverse_gpu[i_f] * pitch);
            p_f[p_step_index_D(i_Global)] = f[i_f]; //- prc(6.) * Density_gpu * w_equ_gpu[i_f] * prc(0.5)*dt_LBM_gpu*body_accelerate_gpu[2] * (T_P)ez_gpu[i_f];
        }
    }
}
//================================================================================================================================================================= =
//----------------------periodic boundary  ----------------------
// for one-step index (OSI) algorithm
//================================================================================================================================================================= =
__global__ void periodic_BC_GPU_face_xy(I_INT imax, T_P *pdf_old_gpu, I_INT pitch, int ntime)
{
    // Indexing (Thread)
    I_INT i = threadIdx.x + blockIdx.x * blockDim.x;
    I_INT ix, iy, iz;
    I_INT i_Global, i_Global_peri; // the index of point in the global distribution function array
    int i_f;                       // column //index of discrete velocities
    T_P *p_f;                      // the point to rowHead(fi in first grid)
    if (i < imax)
    {
        // z direction (top wall)
        ix = i % (NXG0_D) + 1;
        iy = i / (NXG0_D) + 1;
        iz = NZG0_D;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
        i_Global_peri = p_index_D(ix, iy, 0);
// printf("___ %d , %d, %d, %d, %d___", int(BB_xy_top_gpu[0]), int(BB_xy_top_gpu[1]), int(BB_xy_top_gpu[2]), int(BB_xy_top_gpu[3]), int(BB_xy_top_gpu[4]));
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = Reverse_gpu[BB_xy_top_gpu[i3]];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            p_f[p_step_index_D(i_Global)] = p_f[p_step_index_D(i_Global_peri)];
        }
        //__syncthreads();
        // z direction (bottom wall)
        iz = 1;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
        i_Global_peri = p_index_D(ix, iy, NZG0_D + 1);
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = BB_xy_top_gpu[i3];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            p_f[p_step_index_D(i_Global)] = p_f[p_step_index_D(i_Global_peri)];
        }
    }
}
__global__ void periodic_BC_GPU_face_yz(I_INT imax, T_P *pdf_old_gpu, I_INT pitch, int ntime)
{
    // Indexing (Thread)
    I_INT i = threadIdx.x + blockIdx.x * blockDim.x;
    I_INT ix, iy, iz;
    I_INT i_Global, i_Global_peri; // the index of point in the global distribution function array
    int i_f;                       // column //index of discrete velocities
    T_P *p_f;                      // the point to rowHead(fi in first grid)

    if (i < imax)
    {
        // x direction (front wall)
        iy = i % (NYG0_D) + 1;
        iz = i / (NYG0_D) + 1;
        ix = NXG0_D;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
        i_Global_peri = p_index_D(0, iy, iz);
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = Reverse_gpu[BB_yz_front_gpu[i3]];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            p_f[p_step_index_D(i_Global)] = p_f[p_step_index_D(i_Global_peri)];
        }
        // x direction (back wall)
        ix = 1;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
        i_Global_peri = p_index_D(NXG0_D + 1, iy, iz);
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = BB_yz_front_gpu[i3];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            p_f[p_step_index_D(i_Global)] = p_f[p_step_index_D(i_Global_peri)];
        }
    }
}
__global__ void periodic_BC_GPU_face_xz(I_INT imax, T_P *pdf_old_gpu, I_INT pitch, int ntime)
{
    // Indexing (Thread)
    I_INT i = threadIdx.x + blockIdx.x * blockDim.x;
    I_INT ix, iy, iz;
    I_INT i_Global, i_Global_peri; // the index of point in the global distribution function array
    int i_f;                       // column //index of discrete velocities
    T_P *p_f;                      // the point to rowHead(fi in first grid)

    if (i < imax)
    {
        // y direction (Right wall)
        ix = i % (NXG0_D) + 1;
        iz = i / (NXG0_D) + 1;
        iy = NYG0_D;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
        i_Global_peri = p_index_D(ix, 0, iz);
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = Reverse_gpu[BB_xz_Right_gpu[i3]];
            //__syncthreads();
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            p_f[p_step_index_D(i_Global)] = p_f[p_step_index_D(i_Global_peri)];
        }
        // x direction (Left wall)
        iy = 1;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
        i_Global_peri = p_index_D(ix, NYG0_D + 1, iz);
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = BB_xz_Right_gpu[i3];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            p_f[p_step_index_D(i_Global)] = p_f[p_step_index_D(i_Global_peri)];
        }
    }
}
__global__ void periodic_BC_GPU_edge_y(I_INT imax, T_P *pdf_old_gpu, I_INT pitch, int ntime)
{

    int i_f; // column //index of discrete velocities
    // int i_grid ;    //row  //index of grid position coordinates
    T_P *p_f; // the point to rowHead(fi in first grid)

    // Indexing (Thread)
    int i = threadIdx.x + blockIdx.x * blockDim.x;
    int ix, iy, iz;

    // T_P f[NDIR];

    if (i < imax)
    {
        // (1,iy,1)
        iy = i + 1;
        ix = 1;
        iz = 1;
        i_f = 11;
        p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
        p_f[p_step_index_D(p_index_D(ix, iy, iz))] =
            p_f[p_step_index_D(p_index_D(NXG0_D + 1, iy, NZG0_D + 1))];

        // (NXG0_D,iy,1)
        iy = i + 1;
        ix = NXG0_D;
        iz = 1;
        i_f = 12;
        p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
        p_f[p_step_index_D(p_index_D(ix, iy, iz))] =
            p_f[p_step_index_D(p_index_D(0, iy, NZG0_D + 1))];

        // (1,iy,NZG0_D)
        iy = i + 1;
        ix = 1;
        iz = NZG0_D;
        i_f = 13;
        p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
        p_f[p_step_index_D(p_index_D(ix, iy, iz))] =
            p_f[p_step_index_D(p_index_D(NXG0_D + 1, iy, 0))];

        // (NXG0_D,iy,NZG0_D)
        iy = i + 1;
        ix = NXG0_D;
        iz = NZG0_D;
        i_f = 14;
        p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
        p_f[p_step_index_D(p_index_D(ix, iy, iz))] =
            p_f[p_step_index_D(p_index_D(0, iy, 0))];
    }
}
//================================================================================================================================================================= =
//----------------------Zhao-Li et al., 2002, Non-equilibrium extrapolation method for boundary condition for inlet/outlet boundary----------------------
// for one-step index (OSI) algorithm
//================================================================================================================================================================= =
__global__ void BC_Guo2002_GPU_face_xz_save(I_INT imax, T_P *pdf_old_gpu, I_INT pitch, int ntime, T_P *Boundary_xz0_gpu, T_P *Boundary_xz1_gpu)
{            // Zhao-Li et al., 2002, Non-equilibrium extrapolation method
    int i_f; // column //index of discrete velocities
    // int i_grid ;    //row  //index of grid position coordinates
    T_P *p_f; // the point to rowHead(fi in first grid)

    // Indexing (Thread)
    I_INT i = threadIdx.x + blockIdx.x * blockDim.x;
    I_INT ix, iy, iz;
    // Represents the position in the grid that needs to be solved
    // int position;
    if (i < imax)
    {
        // y direction (left : inlet)
        ix = i % (NXG0_D) + 1;
        iy = 1;
        iz = i / (NXG0_D) + 1;
        I_INT i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array

        T_P fbb[NDIR], rhobb, uxbb, uybb, uzbb; // fluid node which neighbour of boundary node
        T_P feqbb, edotubb, udotubb;
        T_P rho = Density_gpu;
        T_P ux = 0, uz = 0, uy = 0.; // boundary node
        T_P feqb, edotu, udotu;

        rhobb = 0;
        uxbb = 0;
        uybb = 0;
        uzbb = 0;
#pragma unroll
        for (int i3 = 0; i3 < NDIR; i3++)
        {
            i_f = i3;
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
#if (GUO_EXTRAPOLATION == BEFORE_COLLSION)
            fbb[i_f] = p_f[p_step_index_D((i_Global))];
#elif (GUO_EXTRAPOLATION == AFTER_COLLSION)
            fbb[i_f] = p_f[p_step_index_D((i_Global) + Length_gpu[i_f])];
#endif
            rhobb = rhobb + fbb[i_f];
            uxbb = uxbb + fbb[i_f] * (T_P)ex_gpu[i_f];
            uybb = uybb + fbb[i_f] * (T_P)ey_gpu[i_f];
            uzbb = uzbb + fbb[i_f] * (T_P)ez_gpu[i_f];
        }
        uxbb = uxbb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[0] * overc_LBM_gpu;
        uybb = uybb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[1] * overc_LBM_gpu;
        uzbb = uzbb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[2] * overc_LBM_gpu;
        udotubb = uxbb * uxbb + uybb * uybb + uzbb * uzbb;

        // rho = Density_gpu * (1. + sqrt(3.)*Velocity_Bound_gpu);//rhobb;
        // ux=uxbb; uz=uybb; uy = uzbb;
        rho = rhobb; // Density_gpu ;//rhobb;
        ux = 0. * overc_LBM_gpu;
        uz = 0. * overc_LBM_gpu;
        uy = Velocity_Bound_gpu; //-Lid_velocity_gpu + (prc(2.)*Lid_velocity_gpu*(iz-ez_gpu[i_f]))/(NZG0_D+1);
        udotu = ux * ux + uy * uy + uz * uz;
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = BB_xz_Right_gpu[i3];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            edotubb = (T_P)ex_gpu[i_f] * uxbb + (T_P)ey_gpu[i_f] * uybb + (T_P)ez_gpu[i_f] * uzbb;
            feqbb = rhobb * w_equ_gpu[i_f] * (prc(1.) + prc(3.) * edotubb + prc(4.5) * edotubb * edotubb - prc(1.5) * udotubb);
            edotu = (T_P)ex_gpu[i_f] * ux + (T_P)ey_gpu[i_f] * uy + (T_P)ez_gpu[i_f] * uz;
            feqb = rho * w_equ_gpu[i_f] * (prc(1.) + prc(3.) * edotu + prc(4.5) * edotu * edotu - prc(1.5) * udotu);
#if (GUO_EXTRAPOLATION == BEFORE_COLLSION)
            Boundary_xz0_gpu[i * 5 + i3] = feqb + (1. - SRT_OverTau_gpu) * (fbb[i_f] - feqbb);
#elif (GUO_EXTRAPOLATION == AFTER_COLLSION)
            p_f[p_step_index_D((i_Global))] = feqb + (fbb[i_f] - feqbb);
#endif
        }

        // y direction (right : outlet)
        iy = NYG0_D;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
        rhobb = 0;
        uxbb = 0;
        uybb = 0;
        uzbb = 0;
#pragma unroll
        for (int i3 = 0; i3 < NDIR; i3++)
        {
            i_f = i3;
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
#if (GUO_EXTRAPOLATION == BEFORE_COLLSION)
            fbb[i_f] = p_f[p_step_index_D((i_Global))];
#elif (GUO_EXTRAPOLATION == AFTER_COLLSION)
            fbb[i_f] = p_f[p_step_index_D((i_Global) + Length_gpu[i_f])];
#endif
            rhobb = rhobb + fbb[i_f];
            uxbb = uxbb + fbb[i_f] * (T_P)ex_gpu[i_f];
            uybb = uybb + fbb[i_f] * (T_P)ey_gpu[i_f];
            uzbb = uzbb + fbb[i_f] * (T_P)ez_gpu[i_f];
        }
        uxbb = uxbb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[0] * overc_LBM_gpu;
        uybb = uybb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[1] * overc_LBM_gpu;
        uzbb = uzbb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[2] * overc_LBM_gpu;
        udotubb = uxbb * uxbb + uybb * uybb + uzbb * uzbb;

        rho = rhobb; // rhobb;//Density_gpu;
        // ux=uxbb; uz=uybb; uy = uzbb;
        // udotu = ux * ux + uy * uy + uz * uz;
        ux = 0. * overc_LBM_gpu;
        uz = 0. * overc_LBM_gpu;
        uy = Velocity_Bound_gpu; //-Lid_velocity_gpu + (prc(2.)*Lid_velocity_gpu*(iz-ez_gpu[i_f]))/(NZG0_D+1);
        udotu = ux * ux + uy * uy + uz * uz;
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = Reverse_gpu[BB_xz_Right_gpu[i3]];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            edotubb = (T_P)ex_gpu[i_f] * uxbb + (T_P)ey_gpu[i_f] * uybb + (T_P)ez_gpu[i_f] * uzbb;
            feqbb = rhobb * w_equ_gpu[i_f] * (prc(1.) + prc(3.) * edotubb + prc(4.5) * edotubb * edotubb - prc(1.5) * udotubb);
            edotu = (T_P)ex_gpu[i_f] * ux + (T_P)ey_gpu[i_f] * uy + (T_P)ez_gpu[i_f] * uz;
            feqb = rho * w_equ_gpu[i_f] * (prc(1.) + prc(3.) * edotu + prc(4.5) * edotu * edotu - prc(1.5) * udotu);
#if (GUO_EXTRAPOLATION == BEFORE_COLLSION)
            Boundary_xz1_gpu[i * 5 + i3] = feqb + (1. - SRT_OverTau_gpu) * (fbb[i_f] - feqbb);
#elif (GUO_EXTRAPOLATION == AFTER_COLLSION)
            p_f[p_step_index_D((i_Global))] = feqb + (fbb[i_f] - feqbb);
#endif
            // feqb + (fbb[i_f] - feqbb);// feqb + (1.-SRT_OverTau_gpu) *(fbb[i_f] - feqbb);
            // printf("\n\t i_f = %f", p_f[p_step_index_D((i_Global) )]);
        }
    }
}
__global__ void BC_Guo2002_GPU_face_xz_load(I_INT imax, T_P *pdf_old_gpu, I_INT pitch, int ntime, T_P *Boundary_xz0_gpu, T_P *Boundary_xz1_gpu)
{ // Zhao-Li et al., 2002, Non-equilibrium extrapolation method

    int i_f; // column //index of discrete velocities
    // int i_grid ;    //row  //index of grid position coordinates
    T_P *p_f; // the point to rowHead(fi in first grid)

    // Indexing (Thread)
    I_INT i = threadIdx.x + blockIdx.x * blockDim.x;
    I_INT ix, iy, iz;
    I_INT i_Global; // the index of point in the global distribution function array
    // Represents the position in the grid that needs to be solved
    // int position;

    if (i < imax)
    {
        ix = i % (NXG0_D) + 1;
        iz = i / (NXG0_D) + 1;
        T_P fw;
        // y direction (left : inlet)
        iy = 1;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = BB_xz_Right_gpu[i3];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            fw = Boundary_xz0_gpu[(i * 5 + i3)];
            p_f[p_step_index_D((i_Global))] = fw;
        }

        // y direction (right : outlet)
        iy = NYG0_D;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = Reverse_gpu[BB_xz_Right_gpu[i3]];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            fw = Boundary_xz1_gpu[(i * 5 + i3)];
            p_f[p_step_index_D((i_Global))] = fw;
        }
    }
}
__global__ void BC_Guo2002_GPU_face_xy_save(I_INT imax, T_P *pdf_old_gpu, I_INT pitch, int ntime, T_P *Boundary_xy0_gpu, T_P *Boundary_xy1_gpu)
{ // Zhao-Li et al., 2002, Non-equilibrium extrapolation method

    int i_f; // column //index of discrete velocities
    // int i_grid ;    //row  //index of grid position coordinates
    T_P *p_f; // the point to rowHead(fi in first grid)

    // Indexing (Thread)
    I_INT i = threadIdx.x + blockIdx.x * blockDim.x;
    I_INT ix, iy, iz;
    I_INT i_Global; // the index of point in the global distribution function array
    // Represents the position in the grid that needs to be solved
    // int position;

    if (i < imax)
    {
        ix = i % (NXG0_D) + 1;
        iy = i / (NXG0_D) + 1;
        iz = 1;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array

        T_P fbb[NDIR], rhobb, uxbb, uybb, uzbb; // fluid node which neighbour of boundary node
        T_P feqbb, edotubb, udotubb;
        T_P rho = Density_gpu;
        T_P ux = 0, uz = 0, uy = 0.; // boundary node
        T_P feqb, edotu, udotu;

        // //z direction (bottom)
        rhobb = 0;
        uxbb = 0;
        uybb = 0;
        uzbb = 0;
#pragma unroll
        for (int i3 = 0; i3 < NDIR; i3++)
        {
            i_f = i3;
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
#if (GUO_EXTRAPOLATION == BEFORE_COLLSION)
            fbb[i_f] = p_f[p_step_index_D((i_Global))];
#elif (GUO_EXTRAPOLATION == AFTER_COLLSION)
            fbb[i_f] = p_f[p_step_index_D((i_Global) + Length_gpu[i_f])];
#endif
            rhobb = rhobb + fbb[i_f];
            uxbb = uxbb + fbb[i_f] * (T_P)ex_gpu[i_f];
            uybb = uybb + fbb[i_f] * (T_P)ey_gpu[i_f];
            uzbb = uzbb + fbb[i_f] * (T_P)ez_gpu[i_f];
        }
        uxbb = uxbb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[0] * overc_LBM_gpu;
        uybb = uybb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[1] * overc_LBM_gpu;
        uzbb = uzbb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[2] * overc_LBM_gpu;
        udotubb = uxbb * uxbb + uybb * uybb + uzbb * uzbb;

        rho = rhobb; // Density_gpu * (1. + sqrt(3.)*Velocity_Bound_gpu);//rhobb;
        ux = 0.;
        uz = 0.;
        uy = Velocity_Bound_gpu;
        udotu = ux * ux + uy * uy + uz * uz;
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = BB_xy_top_gpu[i3];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            edotubb = (T_P)ex_gpu[i_f] * uxbb + (T_P)ey_gpu[i_f] * uybb + (T_P)ez_gpu[i_f] * uzbb;
            feqbb = rhobb * w_equ_gpu[i_f] * (prc(1.) + prc(3.) * edotubb + prc(4.5) * edotubb * edotubb - prc(1.5) * udotubb);
            edotu = (T_P)ex_gpu[i_f] * ux + (T_P)ey_gpu[i_f] * uy + (T_P)ez_gpu[i_f] * uz;
            feqb = rho * w_equ_gpu[i_f] * (prc(1.) + prc(3.) * edotu + prc(4.5) * edotu * edotu - prc(1.5) * udotu);
#if (GUO_EXTRAPOLATION == BEFORE_COLLSION)
            Boundary_xy0_gpu[i * 5 + i3] = feqb + (1. - SRT_OverTau_gpu) * (fbb[i_f] - feqbb);
#elif (GUO_EXTRAPOLATION == AFTER_COLLSION)
            p_f[p_step_index_D((i_Global))] = feqb + (fbb[i_f] - feqbb);
#endif
        }

        // z direction (top)
        iz = NZG0_D;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
        rhobb = 0;
        uxbb = 0;
        uybb = 0;
        uzbb = 0;
#pragma unroll
        for (int i3 = 0; i3 < NDIR; i3++)
        {
            i_f = i3;
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
#if (GUO_EXTRAPOLATION == BEFORE_COLLSION)
            fbb[i_f] = p_f[p_step_index_D((i_Global))];
#elif (GUO_EXTRAPOLATION == AFTER_COLLSION)
            fbb[i_f] = p_f[p_step_index_D((i_Global) + Length_gpu[i_f])];
#endif
            rhobb = rhobb + fbb[i_f];
            uxbb = uxbb + fbb[i_f] * (T_P)ex_gpu[i_f];
            uybb = uybb + fbb[i_f] * (T_P)ey_gpu[i_f];
            uzbb = uzbb + fbb[i_f] * (T_P)ez_gpu[i_f];
        }
        uxbb = uxbb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[0] * overc_LBM_gpu;
        uybb = uybb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[1] * overc_LBM_gpu;
        uzbb = uzbb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[2] * overc_LBM_gpu;
        udotubb = uxbb * uxbb + uybb * uybb + uzbb * uzbb;

        rho = Density_gpu;
        // ux=uxbb; uz=uybb; uy = uzbb;//Lid_velocity_gpu;
        ux = prc(0.);
        uz = prc(0.);
        uy = Velocity_Bound_gpu;
        udotu = ux * ux + uy * uy + uz * uz;
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = Reverse_gpu[BB_xy_top_gpu[i3]];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            edotubb = (T_P)ex_gpu[i_f] * uxbb + (T_P)ey_gpu[i_f] * uybb + (T_P)ez_gpu[i_f] * uzbb;
            feqbb = rhobb * w_equ_gpu[i_f] * (prc(1.) + prc(3.) * edotubb + prc(4.5) * edotubb * edotubb - prc(1.5) * udotubb);
            edotu = (T_P)ex_gpu[i_f] * ux + (T_P)ey_gpu[i_f] * uy + (T_P)ez_gpu[i_f] * uz;
            feqb = rho * w_equ_gpu[i_f] * (prc(1.) + prc(3.) * edotu + prc(4.5) * edotu * edotu - prc(1.5) * udotu);
#if (GUO_EXTRAPOLATION == BEFORE_COLLSION)
            Boundary_xy1_gpu[i * 5 + i3] = feqb + (1. - SRT_OverTau_gpu) * (fbb[i_f] - feqbb);
#elif (GUO_EXTRAPOLATION == AFTER_COLLSION)
            p_f[p_step_index_D((i_Global))] = feqb + (fbb[i_f] - feqbb);
#endif
            // feqb + (fbb[i_f] - feqbb);// feqb + (1.-SRT_OverTau_gpu) *(fbb[i_f] - feqbb);
            // printf("\n\t i_f = %f", p_f[p_step_index_D((i_Global) )]);
        }
    }
}
__global__ void BC_Guo2002_GPU_face_xy_load(I_INT imax, T_P *pdf_old_gpu, I_INT pitch, int ntime, T_P *Boundary_xy0_gpu, T_P *Boundary_xy1_gpu)
{ // Zhao-Li et al., 2002, Non-equilibrium extrapolation method

    int i_f; // column //index of discrete velocities
    // int i_grid ;    //row  //index of grid position coordinates
    T_P *p_f; // the point to rowHead(fi in first grid)

    // Indexing (Thread)
    I_INT i = threadIdx.x + blockIdx.x * blockDim.x;
    I_INT ix, iy, iz;
    I_INT i_Global; // the index of point in the global distribution function array
    // Represents the position in the grid that needs to be solved
    // int position;

    if (i < imax)
    {
        ix = i % (NXG0_D) + 1;
        iy = i / (NXG0_D) + 1;
        T_P fw;
        // z direction (bottom)
        iz = 1;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = BB_xy_top_gpu[i3];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            fw = Boundary_xy0_gpu[(i * 5 + i3)];
            p_f[p_step_index_D((i_Global))] = fw;
        }

        // z direction (top)
        iz = NZG0_D;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = Reverse_gpu[BB_xy_top_gpu[i3]];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            fw = Boundary_xy1_gpu[(i * 5 + i3)];
            p_f[p_step_index_D((i_Global))] = fw;
        }
    }
}
__global__ void BC_Guo2002_GPU_face_yz_save(I_INT imax, T_P *pdf_old_gpu, I_INT pitch, int ntime, T_P *Boundary_yz0_gpu, T_P *Boundary_yz1_gpu)
{ // Zhao-Li et al., 2002, Non-equilibrium extrapolation method

    int i_f; // column //index of discrete velocities
    // int i_grid ;    //row  //index of grid position coordinates
    T_P *p_f; // the point to rowHead(fi in first grid)

    // Indexing (Thread)
    I_INT i = threadIdx.x + blockIdx.x * blockDim.x;
    I_INT ix, iy, iz;
    I_INT i_Global; // the index of point in the global distribution function array
    // Represents the position in the grid that needs to be solved
    // int position;

    if (i < imax)
    {
        // x direction (back)
        ix = 1;
        iy = i % (NYG0_D) + 1;
        iz = i / (NYG0_D) + 1;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array

        T_P fbb[NDIR], rhobb, uxbb, uybb, uzbb; // fluid node which neighbour of boundary node
        T_P feqbb, edotubb, udotubb;
        T_P rho = Density_gpu;
        T_P ux = 0, uz = 0, uy = 0.; // boundary node
        T_P feqb, edotu, udotu;

        rhobb = 0;
        uxbb = 0;
        uybb = 0;
        uzbb = 0;
#pragma unroll
        for (int i3 = 0; i3 < NDIR; i3++)
        {
            i_f = i3;
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
#if (GUO_EXTRAPOLATION == BEFORE_COLLSION)
            fbb[i_f] = p_f[p_step_index_D((i_Global))];
#elif (GUO_EXTRAPOLATION == AFTER_COLLSION)
            fbb[i_f] = p_f[p_step_index_D((i_Global) + Length_gpu[i_f])];
#endif
            rhobb = rhobb + fbb[i_f];
            uxbb = uxbb + fbb[i_f] * (T_P)ex_gpu[i_f];
            uybb = uybb + fbb[i_f] * (T_P)ey_gpu[i_f];
            uzbb = uzbb + fbb[i_f] * (T_P)ez_gpu[i_f];
        }
        uxbb = uxbb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[0] * overc_LBM_gpu;
        uybb = uybb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[1] * overc_LBM_gpu;
        uzbb = uzbb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[2] * overc_LBM_gpu;
        udotubb = uxbb * uxbb + uybb * uybb + uzbb * uzbb;

        // rho = Density_gpu * (1. + sqrt(3.)*Velocity_Bound_gpu);//rhobb;
        // ux=0; uz=0; uy = 0;
        rho = rhobb; // rhobb;
        ux = 0.;
        uz = 0.;
        uy = Velocity_Bound_gpu;
        udotu = ux * ux + uy * uy + uz * uz;
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = BB_yz_front_gpu[i3];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            edotubb = (T_P)ex_gpu[i_f] * uxbb + (T_P)ey_gpu[i_f] * uybb + (T_P)ez_gpu[i_f] * uzbb;
            feqbb = rhobb * w_equ_gpu[i_f] * (prc(1.) + prc(3.) * edotubb + prc(4.5) * edotubb * edotubb - prc(1.5) * udotubb);
            edotu = (T_P)ex_gpu[i_f] * ux + (T_P)ey_gpu[i_f] * uy + (T_P)ez_gpu[i_f] * uz;
            feqb = rho * w_equ_gpu[i_f] * (prc(1.) + prc(3.) * edotu + prc(4.5) * edotu * edotu - prc(1.5) * udotu);
#if (GUO_EXTRAPOLATION == BEFORE_COLLSION)
            Boundary_yz0_gpu[i * 5 + i3] = feqb + (1. - SRT_OverTau_gpu) * (fbb[i_f] - feqbb);
#elif (GUO_EXTRAPOLATION == AFTER_COLLSION)
            p_f[p_step_index_D((i_Global))] = feqb + (fbb[i_f] - feqbb);
#endif
        }

        // x direction (front)
        ix = NXG0_D;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
        rhobb = 0;
        uxbb = 0;
        uybb = 0;
        uzbb = 0;
#pragma unroll
        for (int i3 = 0; i3 < NDIR; i3++)
        {
            i_f = i3;
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
#if (GUO_EXTRAPOLATION == BEFORE_COLLSION)
            fbb[i_f] = p_f[p_step_index_D((i_Global))];
#elif (GUO_EXTRAPOLATION == AFTER_COLLSION)
            fbb[i_f] = p_f[p_step_index_D((i_Global) + Length_gpu[i_f])];
#endif
            rhobb = rhobb + fbb[i_f];
            uxbb = uxbb + fbb[i_f] * (T_P)ex_gpu[i_f];
            uybb = uybb + fbb[i_f] * (T_P)ey_gpu[i_f];
            uzbb = uzbb + fbb[i_f] * (T_P)ez_gpu[i_f];
        }
        uxbb = uxbb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[0] * overc_LBM_gpu;
        uybb = uybb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[1] * overc_LBM_gpu;
        uzbb = uzbb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[2] * overc_LBM_gpu;
        udotubb = uxbb * uxbb + uybb * uybb + uzbb * uzbb;

        rho = rhobb; // rhobb;
        // ux=0; uz=0; uy = -Lid_velocity_gpu+2*Lid_velocity_gpu*(iz)/(NZG0_D+1) ;
        ux = 0.;
        uz = 0.;
        uy = Velocity_Bound_gpu;
        udotu = ux * ux + uy * uy + uz * uz;
// udotu = ux * ux + uy * uy + uz * uz;
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = Reverse_gpu[BB_yz_front_gpu[i3]];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            edotubb = (T_P)ex_gpu[i_f] * uxbb + (T_P)ey_gpu[i_f] * uybb + (T_P)ez_gpu[i_f] * uzbb;
            feqbb = rhobb * w_equ_gpu[i_f] * (prc(1.) + prc(3.) * edotubb + prc(4.5) * edotubb * edotubb - prc(1.5) * udotubb);
            edotu = (T_P)ex_gpu[i_f] * ux + (T_P)ey_gpu[i_f] * uy + (T_P)ez_gpu[i_f] * uz;
            feqb = rho * w_equ_gpu[i_f] * (prc(1.) + prc(3.) * edotu + prc(4.5) * edotu * edotu - prc(1.5) * udotu);
#if (GUO_EXTRAPOLATION == BEFORE_COLLSION)
            Boundary_yz1_gpu[i * 5 + i3] = feqb + (1. - SRT_OverTau_gpu) * (fbb[i_f] - feqbb);
#elif (GUO_EXTRAPOLATION == AFTER_COLLSION)
            p_f[p_step_index_D((i_Global))] = feqb + (fbb[i_f] - feqbb);
#endif
            // feqb + (fbb[i_f] - feqbb);// feqb + (1.-SRT_OverTau_gpu) *(fbb[i_f] - feqbb);
            // printf("\n\t i_f = %f", p_f[p_step_index_D((i_Global) )]);
        }
    }
}
__global__ void BC_Guo2002_GPU_face_yz_load(I_INT imax, T_P *pdf_old_gpu, I_INT pitch, int ntime, T_P *Boundary_yz0_gpu, T_P *Boundary_yz1_gpu)
{ // Zhao-Li et al., 2002, Non-equilibrium extrapolation method

    int i_f; // column //index of discrete velocities
    // int i_grid ;    //row  //index of grid position coordinates
    T_P *p_f; // the point to rowHead(fi in first grid)

    // Indexing (Thread)
    I_INT i = threadIdx.x + blockIdx.x * blockDim.x;
    I_INT ix, iy, iz;
    I_INT i_Global; // the index of point in the global distribution function array
    // Represents the position in the grid that needs to be solved
    // int position;

    if (i < imax)
    {
        iy = i % (NYG0_D) + 1;
        iz = i / (NYG0_D) + 1;
        T_P fw;
        // z direction (back)
        ix = 1;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = BB_yz_front_gpu[i3];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            fw = Boundary_yz0_gpu[(i * 5 + i3)];
            p_f[p_step_index_D((i_Global))] = fw;
        }

        // z direction (front)
        ix = NXG0_D;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = Reverse_gpu[BB_yz_front_gpu[i3]];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            fw = Boundary_yz1_gpu[(i * 5 + i3)];
            p_f[p_step_index_D((i_Global))] = fw;
        }
    }
}
__global__ void BC_NEEM_AfterColl_GPU_face_xz(I_INT imax, T_P *pdf_old_gpu, I_INT pitch, int ntime, T_P *Boundary_xz0_gpu, T_P *Boundary_xz1_gpu)
{             // Non-equilibrium extrapolation method (NEEM) with the disturabtion after collision
    int i_f;  // column //index of discrete velocities
    T_P *p_f; // the point to rowHead(fi in first grid)
    // Indexing (Thread)
    I_INT i = threadIdx.x + blockIdx.x * blockDim.x;
    I_INT ix, iy, iz;
    I_INT i_Global; // the index of point in the global distribution function array
    // Represents the position in the grid that needs to be solved
    // int position;

    if (i < imax)
    {
        // y direction (left : inlet)
        ix = i % (NXG0_D) + 1;
        iy = 1;
        iz = i / (NXG0_D) + 1;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array

        T_P fbb[NDIR], rhobb, uxbb, uybb, uzbb; // fluid node which neighbour of boundary node
        T_P feqbb, edotubb, udotubb;
        T_P rho = Density_gpu;
        T_P ux = 0, uz = 0, uy = Lid_velocity_gpu; // boundary node
        T_P feqb, edotu, udotu;

        rhobb = 0;
        uxbb = 0;
        uybb = 0;
        uzbb = 0;
#pragma unroll
        for (int i3 = 0; i3 < NDIR; i3++)
        {
            i_f = i3;
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            fbb[i_f] = p_f[p_step_index_D((i_Global) + Length_gpu[i_f])];
            rhobb = rhobb + fbb[i_f];
            uxbb = uxbb + fbb[i_f] * (T_P)ex_gpu[i_f];
            uybb = uybb + fbb[i_f] * (T_P)ey_gpu[i_f];
            uzbb = uzbb + fbb[i_f] * (T_P)ez_gpu[i_f];
        }
        uxbb = uxbb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[0] * overc_LBM_gpu;
        uybb = uybb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[1] * overc_LBM_gpu;
        uzbb = uzbb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[2] * overc_LBM_gpu;
        udotubb = uxbb * uxbb + uybb * uybb + uzbb * uzbb;

        // rho = Density_gpu * (1. + sqrt(3.)*Velocity_Bound_gpu);//rhobb;
        rho = rhobb; // Density_gpu ;//rhobb;
        ux = 0;
        uz = 0;
        uy = Velocity_Bound_gpu;
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = BB_xz_Right_gpu[i3];
            ux = 0. * overc_LBM_gpu;
            uz = 0. * overc_LBM_gpu;
            uy = 0. * overc_LBM_gpu; //-Lid_velocity_gpu + (prc(2.)*Lid_velocity_gpu*(iz-ez_gpu[i_f]))/(NZG0_D+1);
            udotu = ux * ux + uy * uy + uz * uz;
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            edotubb = (T_P)ex_gpu[i_f] * uxbb + (T_P)ey_gpu[i_f] * uybb + (T_P)ez_gpu[i_f] * uzbb;
            feqbb = rhobb * w_equ_gpu[i_f] * (prc(1.) + prc(3.) * edotubb + prc(4.5) * edotubb * edotubb - prc(1.5) * udotubb);
            edotu = (T_P)ex_gpu[i_f] * ux + (T_P)ey_gpu[i_f] * uy + (T_P)ez_gpu[i_f] * uz;
            feqb = rho * w_equ_gpu[i_f] * (prc(1.) + prc(3.) * edotu + prc(4.5) * edotu * edotu - prc(1.5) * udotu);
            p_f[p_step_index_D((i_Global))] = feqb + (fbb[i_f] - feqbb);
        }

        // y direction (right : outlet)
        iy = NYG0_D;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
        rhobb = 0;
        uxbb = 0;
        uybb = 0;
        uzbb = 0;
#pragma unroll
        for (int i3 = 0; i3 < NDIR; i3++)
        {
            i_f = i3;
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            fbb[i_f] = p_f[p_step_index_D((i_Global) + Length_gpu[i_f])];
            rhobb = rhobb + fbb[i_f];
            uxbb = uxbb + fbb[i_f] * (T_P)ex_gpu[i_f];
            uybb = uybb + fbb[i_f] * (T_P)ey_gpu[i_f];
            uzbb = uzbb + fbb[i_f] * (T_P)ez_gpu[i_f];
        }
        uxbb = uxbb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[0] * overc_LBM_gpu;
        uybb = uybb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[1] * overc_LBM_gpu;
        uzbb = uzbb / rhobb + prc(0.5) * dt_LBM_gpu * body_accelerate_gpu[2] * overc_LBM_gpu;
        udotubb = uxbb * uxbb + uybb * uybb + uzbb * uzbb;

        rho = rhobb; // rhobb;//Density_gpu;
// ux=uxbb; uz=uybb; uy = uzbb;
// udotu = ux * ux + uy * uy + uz * uz;
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = Reverse_gpu[BB_xz_Right_gpu[i3]];
            ux = 0. * overc_LBM_gpu;
            uz = 0. * overc_LBM_gpu;
            uy = 0. * overc_LBM_gpu; //-Lid_velocity_gpu + (prc(2.)*Lid_velocity_gpu*(iz-ez_gpu[i_f]))/(NZG0_D+1);
            udotu = ux * ux + uy * uy + uz * uz;
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            edotubb = (T_P)ex_gpu[i_f] * uxbb + (T_P)ey_gpu[i_f] * uybb + (T_P)ez_gpu[i_f] * uzbb;
            feqbb = rhobb * w_equ_gpu[i_f] * (prc(1.) + prc(3.) * edotubb + prc(4.5) * edotubb * edotubb - prc(1.5) * udotubb);
            edotu = (T_P)ex_gpu[i_f] * ux + (T_P)ey_gpu[i_f] * uy + (T_P)ez_gpu[i_f] * uz;
            feqb = rho * w_equ_gpu[i_f] * (prc(1.) + prc(3.) * edotu + prc(4.5) * edotu * edotu - prc(1.5) * udotu);
            p_f[p_step_index_D((i_Global))] = feqb + (fbb[i_f] - feqbb);
        }
    }
}
__global__ void Out_ABC_xz(I_INT imax, T_P *pdf_old_gpu, I_INT pitch, int ntime, T_P *Boundary_xz0_gpu, T_P *Boundary_xz1_gpu)
{ // Zhao-Li et al., 2002, Non-equilibrium extrapolation method

    int i_f; // column //index of discrete velocities
    // int i_grid ;    //row  //index of grid position coordinates
    T_P *p_f; // the point to rowHead(fi in first grid)

    // Indexing (Thread)
    I_INT i = threadIdx.x + blockIdx.x * blockDim.x;
    I_INT ix, iy, iz;
    I_INT i_Global; // the index of point in the global distribution function array
    // Represents the position in the grid that needs to be solved
    // int position;

    if (i < imax)
    {
        // y direction (right : outlet)
        ix = i % (NXG0_D) + 1;
        iz = i / (NXG0_D) + 1;
        T_P fw;
        iy = NYG0_D;
        i_Global = p_index_D(ix, iy, iz); // the index of point in the global distribution function array
#pragma unroll
        for (int i3 = 0; i3 < 5; i3++)
        {
            i_f = Reverse_gpu[BB_xz_Right_gpu[i3]];
            p_f = (T_P *)((char *)pdf_old_gpu + i_f * pitch);
            // fw = Density_gpu * w_equ_gpu[i_f] * (prc(1.) + prc(3.)*(T_P)ey_gpu[i_f]*Lid_velocity_gpu + prc(4.5)*(T_P)ey_gpu[i_f]*Lid_velocity_gpu*(T_P)ey_gpu[i_f]*Lid_velocity_gpu - prc(1.5)* Lid_velocity_gpu*Lid_velocity_gpu);
            fw = Density_gpu * w_equ_gpu[i_f] * (prc(1.));
            p_f[p_step_index_D((i_Global))] = fw;
        }
    }
}
